ATAC-seq accessibility exploration

library(ArchR)

## Loading required package: ggplot2

## Loading required package: SummarizedExperiment

## Loading required package: GenomicRanges

## Loading required package: stats4

## Loading required package: BiocGenerics

## Loading required package: parallel

## 
## Attaching package: 'BiocGenerics'

## The following objects are masked from 'package:parallel':
## 
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, parApply, parCapply, parLapply,
##     parLapplyLB, parRapply, parSapply, parSapplyLB

## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs

## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, basename, cbind, colnames,
##     dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
##     grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
##     order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
##     rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
##     union, unique, unsplit, which, which.max, which.min

## Loading required package: S4Vectors

## 
## Attaching package: 'S4Vectors'

## The following object is masked from 'package:base':
## 
##     expand.grid

## Loading required package: IRanges

## Loading required package: GenomeInfoDb

## Loading required package: Biobase

## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.

## Loading required package: DelayedArray

## Loading required package: matrixStats

## 
## Attaching package: 'matrixStats'

## The following objects are masked from 'package:Biobase':
## 
##     anyMissing, rowMedians

## Loading required package: BiocParallel

## 
## Attaching package: 'DelayedArray'

## The following objects are masked from 'package:matrixStats':
## 
##     colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges

## The following objects are masked from 'package:base':
## 
##     aperm, apply, rowsum

## Loading required package: data.table

## 
## Attaching package: 'data.table'

## The following object is masked from 'package:SummarizedExperiment':
## 
##     shift

## The following object is masked from 'package:GenomicRanges':
## 
##     shift

## The following object is masked from 'package:IRanges':
## 
##     shift

## The following objects are masked from 'package:S4Vectors':
## 
##     first, second

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

## The following object is masked from 'package:S4Vectors':
## 
##     expand

## Loading required package: rhdf5

## Loading required package: magrittr

archProj <- loadArchRProject("../data/scATAC/archR_subHybridClusters/",
                 showLogo = FALSE)

## Successfully loaded ArchRProject!

plotEmbedding(ArchRProj = archProj, colorBy = "cellColData", 
              name = "predAnn", embedding = "UMAP_cor")

## ArchR logging to : ArchRLogs/ArchR-plotEmbedding-c6f51b0612db-Date-2021-12-25_Time-17-20-39.log
## If there is an issue, please report to github with logFile!

## Getting UMAP Embedding

## ColorBy = cellColData

## Plotting Embedding

## 1 
## ArchR logging successful to : ArchRLogs/ArchR-plotEmbedding-c6f51b0612db-Date-2021-12-25_Time-17-20-39.log

archProj <- addMotifAnnotations(ArchRProj = archProj, 
                                motifSet = "cisbp", 
                                name = "Motif",
                                force=TRUE)

## No methods found in package 'IRanges' for request: 'score' when loading 'TFBSTools'
## ArchR logging to : ArchRLogs/ArchR-addMotifAnnotations-c6f52bbff5b9-Date-2021-12-25_Time-17-20-42.log
## If there is an issue, please report to github with logFile!
## peakAnnotation name already exists! Overriding.
## 2021-12-25 17:20:47 : Gettting Motif Set, Species : Mus musculus, 0.001 mins elapsed.
## Using version 2 motifs!
## 2021-12-25 17:20:51 : Finding Motif Positions with motifmatchr!, 0.079 mins elapsed.
## 2021-12-25 17:24:13 : Creating Motif Overlap Matrix, 3.436 mins elapsed.
## 2021-12-25 17:24:16 : Finished Getting Motif Info!, 3.499 mins elapsed.
## ArchR logging successful to : ArchRLogs/ArchR-addMotifAnnotations-c6f52bbff5b9-Date-2021-12-25_Time-17-20-42.log

umapEmbedding <- getEmbedding(archProj, embedding="UMAP_cor")
colnames(umapEmbedding) <- paste0("UMAP",1:2)
umapEmbedding$cellState <- archProj$predAnn
umapEmbedding$cellState[umapEmbedding$cellState == "C1_Inj"] <- "injured"
umapEmbedding$cellState[umapEmbedding$cellState == "C2_Hybrid"] <- "Hybrid"
umapEmbedding$cellState <- factor(umapEmbedding$cellState)
ggplot(umapEmbedding, aes(x=UMAP1, y=UMAP2, col=cellState)) +
  geom_point(size=1/2) +
  theme_classic() 

ggsave("../plots/hbcCellState.png")

## Saving 7 x 5 in image

## get ArchR peak markers for motif enrichment
markerPeaks_subHybrid <- getMarkerFeatures(
    ArchRProj = archProj, 
    useMatrix = "PeakMatrix", 
    groupBy = "predAnn",
    bias = c("TSSEnrichment", "log10(nFrags)"),
    testMethod = "wilcoxon")

## ArchR logging to : ArchRLogs/ArchR-getMarkerFeatures-c6f5479d5b1f-Date-2021-12-25_Time-17-24-21.log
## If there is an issue, please report to github with logFile!
## MatrixClass = Sparse.Integer.Matrix
## 2021-12-25 17:24:22 : Matching Known Biases, 0.004 mins elapsed.
## 2021-12-25 17:24:23 : Computing Pairwise Tests (1 of 5), 0.035 mins elapsed.
## 2021-12-25 17:24:47 : Computing Pairwise Tests (2 of 5), 0.429 mins elapsed.
## 2021-12-25 17:25:09 : Computing Pairwise Tests (3 of 5), 0.791 mins elapsed.
## 2021-12-25 17:25:28 : Computing Pairwise Tests (4 of 5), 1.119 mins elapsed.
## 2021-12-25 17:25:49 : Computing Pairwise Tests (5 of 5), 1.456 mins elapsed.
## ###########
## 2021-12-25 17:26:11 : Completed Pairwise Tests, 1.825 mins elapsed.
## ###########
## ArchR logging successful to : ArchRLogs/ArchR-getMarkerFeatures-c6f5479d5b1f-Date-2021-12-25_Time-17-24-21.log

peakMarkers_subHybrid <- getMarkers(markerPeaks_subHybrid, 
                                    cutOff = "FDR <= 0.01 & Log2FC >= 1.25", 
                                    returnGR = TRUE)
heatmapMarkerPeaks <- plotMarkerHeatmap(
  seMarker = markerPeaks_subHybrid, 
  cutOff = "FDR <= 0.01 & Log2FC >= 1.25", 
  transpose = TRUE
)

## ArchR logging to : ArchRLogs/ArchR-plotMarkerHeatmap-c6f52bc974f4-Date-2021-12-25_Time-17-26-13.log
## If there is an issue, please report to github with logFile!
## Identified 20694 markers!
## Preparing Main Heatmap..
## ArchR logging successful to : ArchRLogs/ArchR-plotMarkerHeatmap-c6f52bc974f4-Date-2021-12-25_Time-17-26-13.log

ComplexHeatmap::draw(heatmapMarkerPeaks, heatmap_legend_side = "bot", annotation_legend_side = "bot")

Check accessibility of Egfr target genes

egfrTargetDf <- read.csv("../data/genesTargets/egfr_targets_unfiltered.csv")
pCluster <- plotEmbedding(ArchRProj = archProj, 
                    colorBy = "cellColData", 
                    name = "predAnn", 
                    embedding = "UMAP_cor")

## ArchR logging to : ArchRLogs/ArchR-plotEmbedding-c6f533d9ffb1-Date-2021-12-25_Time-17-26-20.log
## If there is an issue, please report to github with logFile!

## Getting UMAP Embedding

## ColorBy = cellColData

## Plotting Embedding

## 1 
## ArchR logging successful to : ArchRLogs/ArchR-plotEmbedding-c6f533d9ffb1-Date-2021-12-25_Time-17-26-20.log

pCluster

allGenes <- unname(unique(mcols(archProj@geneAnnotation$genes)$symbol))
egfrTargets <- intersect(as.character(egfrTargetDf$target), allGenes)

## get peak count matrix
ArrowFiles <- getArrowFiles(archProj)
featureDF <- ArchR:::.getFeatureDF(head(ArrowFiles, 2), "PeakMatrix")
peakMat <- ArchR:::.getPartialMatrix(ArrowFiles, 
            featureDF = featureDF, threads = 1, useMatrix = "PeakMatrix", 
            cellNames = rownames(archProj@cellColData), progress = FALSE)

## 2021-12-25 17:26:21 : Getting Partial Matrix 1 of 6, 0 mins elapsed.
## 2021-12-25 17:26:30 : Getting Partial Matrix 2 of 6, 0.156 mins elapsed.
## 2021-12-25 17:26:35 : Getting Partial Matrix 3 of 6, 0.23 mins elapsed.
## 2021-12-25 17:26:39 : Getting Partial Matrix 4 of 6, 0.298 mins elapsed.
## 2021-12-25 17:26:42 : Getting Partial Matrix 5 of 6, 0.356 mins elapsed.
## 2021-12-25 17:26:46 : Getting Partial Matrix 6 of 6, 0.41 mins elapsed.
## 2021-12-25 17:26:54 : Successfully Created Partial Matrix, 0.557 mins elapsed.

## for each gene, get the peak closest to TSS
peakGR <- archProj@peakSet
tssPeaks <- c()
for(gg in 1:length(egfrTargets)){
  curGene <- egfrTargets[gg]
  curID <- which(mcols(peakGR)$nearestGene == curGene)
  if(length(curID) > 0 ){
    curTSSPeak <- curID[which.min(mcols(peakGR)$distToTSS[curID])]
    tssPeaks[gg] <- curTSSPeak
  } else {
    tssPeaks[gg] <- NA
    next
  }
}

egfrTargets <- egfrTargets[!is.na(tssPeaks)]
tssPeaks <- tssPeaks[!is.na(tssPeaks)]
egfrTargetPeaks <- peakMat[tssPeaks,]
# egfrTargetPeakScaledSum <- colSums(t(scale(t(egfrTargetPeaks))))
egfrTargetPeakScaledSum <- colSums(egfrTargetPeaks)


umapEmbedding <- getEmbedding(archProj, embedding="UMAP_cor")
colnames(umapEmbedding) <- paste0("UMAP",1:2)
umapEmbedding$egfrTargetScaledSum <- egfrTargetPeakScaledSum

pEgfrTargets <- ggplot(umapEmbedding, aes(x=UMAP1,
                          y=UMAP2,
                          col=egfrTargetScaledSum)) +
  theme_classic() +
  geom_point() +
  scale_color_gradientn(colors=wesanderson::wes_palette("Zissou1", n=100, type="continuous"))

pEgfrTargets

Check accessibility of peaks containining Ets TFBS motif

motifMatches <- readRDS(archProj@peakAnnotation$Motif$Matches)
etsMotifId <- grep(x=colnames(motifMatches), pattern="Ets")
etsMotifs <- grep(x=colnames(motifMatches), pattern="Ets", value=TRUE)
peaksWithEtsMotifs <- which(rowSums(assays(motifMatches[,etsMotifs])$matches) > 0)
length(peaksWithEtsMotifs)

## [1] 10594

umapEmbedding$etsMotifPeaks <- colSums(peakMat[peaksWithEtsMotifs,])

pEtsMotifs <- ggplot(umapEmbedding, aes(x=UMAP1,
                          y=UMAP2,
                          col=etsMotifPeaks)) +
  theme_classic() +
  geom_point() +
  scale_color_gradientn(colors=wesanderson::wes_palette("Zissou1", n=100, type="continuous"))

pEtsMotifs

Association between Egfr targets’ accessibility and Ets motif peaks’ accessibility

plot(umapEmbedding$egfrTargetScaledSum, umapEmbedding$etsMotifPeaks,
     xlab = "Egfr targets' accessibility",
     ylab = "Ets-motif-containing peaks' accessibility")

Gene markers for (sub)hybrid clusters

plotEmbedding(ArchRProj = archProj, colorBy = "cellColData", 
              name = "predAnn", embedding = "UMAP_cor")

## ArchR logging to : ArchRLogs/ArchR-plotEmbedding-c6f533ad5e17-Date-2021-12-25_Time-17-26-59.log
## If there is an issue, please report to github with logFile!

## Getting UMAP Embedding

## ColorBy = cellColData

## Plotting Embedding

## 1 
## ArchR logging successful to : ArchRLogs/ArchR-plotEmbedding-c6f533ad5e17-Date-2021-12-25_Time-17-26-59.log

## get ArchR gene markers for motif enrichment
markerGenes_subHybrid <- getMarkerFeatures(
    ArchRProj = archProj, 
    useMatrix = "GeneScoreMatrix", 
    groupBy = "predAnn",
    bias = c("TSSEnrichment", "log10(nFrags)"),
    testMethod = "wilcoxon")

## ArchR logging to : ArchRLogs/ArchR-getMarkerFeatures-c6f5b8cf30f-Date-2021-12-25_Time-17-27-01.log
## If there is an issue, please report to github with logFile!

## MatrixClass = Sparse.Double.Matrix

## 2021-12-25 17:27:01 : Matching Known Biases, 0.001 mins elapsed.

## 2021-12-25 17:27:02 : Computing Pairwise Tests (1 of 5), 0.024 mins elapsed.

## 2021-12-25 17:27:24 : Computing Pairwise Tests (2 of 5), 0.386 mins elapsed.

## 2021-12-25 17:27:44 : Computing Pairwise Tests (3 of 5), 0.727 mins elapsed.

## 2021-12-25 17:28:05 : Computing Pairwise Tests (4 of 5), 1.064 mins elapsed.

## 2021-12-25 17:28:23 : Computing Pairwise Tests (5 of 5), 1.376 mins elapsed.

## ###########
## 2021-12-25 17:28:43 : Completed Pairwise Tests, 1.71 mins elapsed.
## ###########

## ArchR logging successful to : ArchRLogs/ArchR-getMarkerFeatures-c6f5b8cf30f-Date-2021-12-25_Time-17-27-01.log

geneMarkers_subHybrid <- getMarkers(markerGenes_subHybrid, 
                                    cutOff = "FDR <= 0.01 & Log2FC >= 1.25")

write.csv(geneMarkers_subHybrid$C2_Hybrid, file="../data/scATAC/markerGenesHybrid.csv")
write.csv(geneMarkers_subHybrid$hybrid1, file="../data/scATAC/markerGenesHybrid1.csv")
write.csv(geneMarkers_subHybrid$hybrid2, file="../data/scATAC/markerGenesHybrid2.csv")
write.csv(geneMarkers_subHybrid$C1_Inj, file="../data/scATAC/markerGenesInjured.csv")
write.csv(geneMarkers_subHybrid$resting, file="../data/scATAC/markerGenesResting.csv")

Marker genes for hybrid1

markerGenesHybrid1 <- geneMarkers_subHybrid$hybrid1$name[1:18]
p <- plotEmbedding(
    ArchRProj = archProj, 
    colorBy = "GeneScoreMatrix", 
    name = markerGenesHybrid1, 
    embedding = "UMAP_cor",
    quantCut = c(0.01, 0.95),
    size=1/3,
    plotAs="points"
)

## Getting ImputeWeights

## No imputeWeights found, returning NULL

## ArchR logging to : ArchRLogs/ArchR-plotEmbedding-c6f56ed12a78-Date-2021-12-25_Time-17-28-44.log
## If there is an issue, please report to github with logFile!

## Getting UMAP Embedding

## ColorBy = GeneScoreMatrix

## Getting Matrix Values...

## 2021-12-25 17:28:44 :

## 1 2 3 4 5 6 
## Plotting Embedding
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 
## ArchR logging successful to : ArchRLogs/ArchR-plotEmbedding-c6f56ed12a78-Date-2021-12-25_Time-17-28-44.log

# clean
p2Hybrid1 <- lapply(p, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})
cowplot::plot_grid(plotlist=p2Hybrid1[1:9])

cowplot::plot_grid(plotlist=p2Hybrid1[10:18])

p2Hybrid1

## $Epas1

## 
## $Ptgfr

## 
## $Kcnj1

## 
## $Tmc5

## 
## $Pla2g2c

## 
## $Muc5b

## 
## $B3galt1

## 
## $Tacr1

## 
## $Sult1c1

## 
## $Arhgef38

## 
## $Cyp4a12a

## 
## $Zdhhc19

## 
## $Spint3

## 
## $Acnat1

## 
## $`9630013K17Rik`

## 
## $Car8

## 
## $Sftpb

## 
## $Tff2

## make plots manually
umapEmbedding <- getEmbedding(archProj, embedding="UMAP_cor")
colnames(umapEmbedding) <- paste0("UMAP",1:2)
featureDF <- ArchR:::.getFeatureDF(head(ArrowFiles, 2), "GeneScoreMatrix")
featureDFHybrid1 <- featureDF[which(featureDF$name %in% markerGenesHybrid1),]
geneActMatHybrid1 <- ArchR:::.getPartialMatrix(ArrowFiles, 
            featureDF = featureDFHybrid1, threads = 1, useMatrix = "GeneScoreMatrix", 
            cellNames = rownames(archProj@cellColData), progress = FALSE)

## 2021-12-25 17:29:11 : Getting Partial Matrix 1 of 6, 0 mins elapsed.
## 2021-12-25 17:29:12 : Getting Partial Matrix 2 of 6, 0.01 mins elapsed.
## 2021-12-25 17:29:12 : Getting Partial Matrix 3 of 6, 0.02 mins elapsed.
## 2021-12-25 17:29:13 : Getting Partial Matrix 4 of 6, 0.026 mins elapsed.
## 2021-12-25 17:29:13 : Getting Partial Matrix 5 of 6, 0.033 mins elapsed.
## 2021-12-25 17:29:13 : Getting Partial Matrix 6 of 6, 0.039 mins elapsed.
## 2021-12-25 17:29:15 : Successfully Created Partial Matrix, 0.066 mins elapsed.

customAlpha <- function(y){
  alpha1 <- pmin(as.numeric(factor(dayjob::colby(y, g=3)))-.8, 1)
  return(alpha1)
}

plistIkHybrid1 <- list()
for(gg in 1:length(markerGenesHybrid1)){
  curDF <- umapEmbedding
  curDF$markerGene <- log1p(geneActMatHybrid1[gg,])
  curP <- ggplot(curDF, aes(x=UMAP1,
                            y=UMAP2,
                            col=markerGene)) +
    theme_classic() +
    ggtitle(markerGenesHybrid1[gg]) +
    geom_point(alpha=pmin(1,(log1p(curDF$markerGene) / max(log1p(curDF$markerGene))) + 0.05)
               , pch=16) +
    scale_color_gradientn(colors=wesanderson::wes_palette("Zissou1", n=100, type="continuous"), breaks=c(0, quantile(curDF$markerGene, probs=seq(0.05, 0.95, length=98)) , max(curDF$markerGene))) +
    theme(legend.position = "none")
  plistIkHybrid1[[gg]] <- curP
}
cowplot::plot_grid(plotlist=plistIkHybrid1[1:9])

cowplot::plot_grid(plotlist=plistIkHybrid1[10:18])

Injured/uninjured contribution of marker genes

geneScoreMat <- ArchR::getMatrixFromProject(archProj,
                            useMatrix = "GeneScoreMatrix")

## ArchR logging to : ArchRLogs/ArchR-getMatrixFromProject-c6f5456bc901-Date-2021-12-25_Time-17-29-31.log
## If there is an issue, please report to github with logFile!

## 2021-12-25 17:29:50 : Organizing colData, 0.317 mins elapsed.

## 2021-12-25 17:29:51 : Organizing rowData, 0.318 mins elapsed.

## 2021-12-25 17:29:51 : Organizing Assays (1 of 1), 0.318 mins elapsed.

## 2021-12-25 17:29:52 : Constructing SummarizedExperiment, 0.35 mins elapsed.

## 2021-12-25 17:29:53 : Finished Matrix Creation, 0.367 mins elapsed.

geneScoreMatHybrid1 <- geneScoreMat[,colData(geneScoreMat)$predAnn == "hybrid1"]

table(geneScoreMatHybrid1$treat)

## 
## Inj  UI 
##  77 120

rafalib::mypar(mfrow=c(3,3))
for(gg in 1:length(markerGenesHybrid1)){
  boxplot(log1p(assays(geneScoreMatHybrid1)$GeneScoreMatrix[which(rowData(geneScoreMatHybrid1)$name == markerGenesHybrid1[gg]),]) ~ geneScoreMatHybrid1$treat,
        xlab = "Treatment", ylab = "Log Gene activity score + 1",
        main = markerGenesHybrid1[gg])
}

Marker genes for hybrid2

markerGenesHybrid2 <- geneMarkers_subHybrid$hybrid2$name[1:18]
p <- plotEmbedding(
    ArchRProj = archProj, 
    colorBy = "GeneScoreMatrix", 
    name = markerGenesHybrid2, 
    embedding = "UMAP_cor",
    quantCut = c(0.01, 0.95),
    size=1/3,
    plotAs="points"
)

## Getting ImputeWeights

## No imputeWeights found, returning NULL

## ArchR logging to : ArchRLogs/ArchR-plotEmbedding-c6f565fa1d3c-Date-2021-12-25_Time-17-29-54.log
## If there is an issue, please report to github with logFile!

## Getting UMAP Embedding

## ColorBy = GeneScoreMatrix

## Getting Matrix Values...

## 2021-12-25 17:29:55 :

## 1 2 3 4 5 6 
## Plotting Embedding
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 
## ArchR logging successful to : ArchRLogs/ArchR-plotEmbedding-c6f565fa1d3c-Date-2021-12-25_Time-17-29-54.log

# clean
p2Hybrid2 <- lapply(p, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})
cowplot::plot_grid(plotlist=p2Hybrid2[1:9])

cowplot::plot_grid(plotlist=p2Hybrid2[10:18])

p2Hybrid2

## $Gm32865

## 
## $`2610016A17Rik`

## 
## $Pax9

## 
## $Pax3

## 
## $Atp10a

## 
## $Kcnt2

## 
## $Aldh1a7

## 
## $Aldh1a1

## 
## $Cfh

## 
## $Cap2

## 
## $Pcare

## 
## $Slc25a21

## 
## $Vwc2

## 
## $Irx5

## 
## $C730002L08Rik

## 
## $Crnde

## 
## $Tmem45a2

## 
## $Psg22

## make plots manually
umapEmbedding <- getEmbedding(archProj, embedding="UMAP_cor")
colnames(umapEmbedding) <- paste0("UMAP",1:2)
featureDF <- ArchR:::.getFeatureDF(head(ArrowFiles, 2), "GeneScoreMatrix")
featureDFHybrid2 <- featureDF[which(featureDF$name %in% markerGenesHybrid2),]
geneActMatHybrid2 <- ArchR:::.getPartialMatrix(ArrowFiles, 
            featureDF = featureDFHybrid2, threads = 1, useMatrix = "GeneScoreMatrix", 
            cellNames = rownames(archProj@cellColData), progress = FALSE)

## 2021-12-25 17:30:17 : Getting Partial Matrix 1 of 6, 0 mins elapsed.
## 2021-12-25 17:30:18 : Getting Partial Matrix 2 of 6, 0.01 mins elapsed.
## 2021-12-25 17:30:18 : Getting Partial Matrix 3 of 6, 0.016 mins elapsed.
## 2021-12-25 17:30:19 : Getting Partial Matrix 4 of 6, 0.024 mins elapsed.
## 2021-12-25 17:30:19 : Getting Partial Matrix 5 of 6, 0.029 mins elapsed.
## 2021-12-25 17:30:19 : Getting Partial Matrix 6 of 6, 0.033 mins elapsed.
## 2021-12-25 17:30:21 : Successfully Created Partial Matrix, 0.068 mins elapsed.

customAlpha <- function(y){
  alpha1 <- pmin(as.numeric(factor(dayjob::colby(y, g=3)))-.8, 1)
  return(alpha1)
}

plistIkHybrid2 <- list()
for(gg in 1:length(markerGenesHybrid2)){
  curDF <- umapEmbedding
  curDF$markerGene <- log1p(geneActMatHybrid2[gg,])
  curP <- ggplot(curDF, aes(x=UMAP1,
                            y=UMAP2,
                            col=markerGene)) +
    theme_classic() +
    ggtitle(markerGenesHybrid2[gg]) +
    geom_point(alpha=pmin(1,(log1p(curDF$markerGene) / max(log1p(curDF$markerGene))) + 0.05)
               , pch=16) +
    scale_color_gradientn(colors=wesanderson::wes_palette("Zissou1", n=100, type="continuous"), breaks=c(0, quantile(curDF$markerGene, probs=seq(0.05, 0.95, length=98)) , max(curDF$markerGene))) +
    theme(legend.position = "none")
  plistIkHybrid2[[gg]] <- curP
}
cowplot::plot_grid(plotlist=plistIkHybrid2[1:9])

cowplot::plot_grid(plotlist=plistIkHybrid2[10:18])

Injured/uninjured contribution of marker genes

geneScoreMatHybrid2 <- geneScoreMat[,colData(geneScoreMat)$predAnn == "hybrid2"]

table(geneScoreMatHybrid2$treat)

## 
## Inj  UI 
##  47 176

rafalib::mypar(mfrow=c(3,3))
for(gg in 1:length(markerGenesHybrid2)){
  boxplot(log1p(assays(geneScoreMatHybrid2)$GeneScoreMatrix[which(rowData(geneScoreMatHybrid2)$name == markerGenesHybrid2[gg]),]) ~ geneScoreMatHybrid2$treat,
        xlab = "Treatment", ylab = "Log Gene activity score + 1",
        main = markerGenesHybrid2[gg])
}

Marker genes for other hybrid cells

markerGenesHybrid <- geneMarkers_subHybrid$C2_Hybrid$name[1:18]
p <- plotEmbedding(
    ArchRProj = archProj, 
    colorBy = "GeneScoreMatrix", 
    name = markerGenesHybrid, 
    embedding = "UMAP_cor",
    quantCut = c(0.01, 0.95),
    size=1/3,
    plotAs="points"
)

## Getting ImputeWeights

## No imputeWeights found, returning NULL

## ArchR logging to : ArchRLogs/ArchR-plotEmbedding-c6f52f990314-Date-2021-12-25_Time-17-30-37.log
## If there is an issue, please report to github with logFile!

## Getting UMAP Embedding

## ColorBy = GeneScoreMatrix

## Getting Matrix Values...

## 2021-12-25 17:30:37 :

## 1 2 3 4 5 6 
## Plotting Embedding
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 
## ArchR logging successful to : ArchRLogs/ArchR-plotEmbedding-c6f52f990314-Date-2021-12-25_Time-17-30-37.log

# clean
p2Hybrid <- lapply(p, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})
cowplot::plot_grid(plotlist=p2Hybrid[1:9])

cowplot::plot_grid(plotlist=p2Hybrid[10:18])

Accesssibility of key lineage-specific genes from other cell types / lineages

# 
markerGenesLineageSpecific <- c("Sox2", "Ascl1", "Kit", "Neurog1", "Neurod1")
p <- plotEmbedding(
    ArchRProj = archProj, 
    colorBy = "GeneScoreMatrix", 
    name = markerGenesLineageSpecific, 
    embedding = "UMAP_cor",
    quantCut = c(0.01, 0.95),
    size=1/3,
    plotAs="points"
)

## Getting ImputeWeights

## No imputeWeights found, returning NULL

## ArchR logging to : ArchRLogs/ArchR-plotEmbedding-c6f552b937f2-Date-2021-12-25_Time-17-30-54.log
## If there is an issue, please report to github with logFile!

## Getting UMAP Embedding

## ColorBy = GeneScoreMatrix

## Getting Matrix Values...

## 2021-12-25 17:30:55 :

## 1 2 3 4 5 6 
## Plotting Embedding
## 1 2 3 4 5 
## ArchR logging successful to : ArchRLogs/ArchR-plotEmbedding-c6f552b937f2-Date-2021-12-25_Time-17-30-54.log

# clean
p2LineageSpecific <- lapply(p, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})
cowplot::plot_grid(plotlist=p2LineageSpecific)

Accessibility of hybrid markes identified using RNA-seq

# marker genes identified in lineage-traced scRNA-seq when combining hybrid1/hybrid2
markerHybridAll <- c("Ly6a", "Adh7", "Mgp", "Selenbp1", "Sftpd", "Sult1d1",
                     "Cyp2f2", "Reg3g", "Krt15", "Bpifa1")
markerHyrbid1 <- c("Sult1d1", "Serpina3n", "Pglyrp1", "Adh7", "Gsta3", "Prdx6",
                   "Reg3g", "Lcn2", "Bpifa1", "Ltf")
markerHyrbid2 <- c("Adh1", "Mgp", "Ly6a", "Ly6d", "Sftpd", "Selenbp1",
                   "Krt15", "Aldh3a1", "Cyp2f2")

Markers for all hybrid cells

pMarkerHybridAll <- plotEmbedding(
    ArchRProj = archProj, 
    colorBy = "GeneScoreMatrix", 
    name = markerHybridAll[1:9], 
    embedding = "UMAP_cor",
    quantCut = c(0.01, 0.95),
    size=1/3,
    plotAs="points"
)

## Getting ImputeWeights

## No imputeWeights found, returning NULL

## ArchR logging to : ArchRLogs/ArchR-plotEmbedding-c6f5d16a29b-Date-2021-12-25_Time-17-31-01.log
## If there is an issue, please report to github with logFile!

## Getting UMAP Embedding

## ColorBy = GeneScoreMatrix

## Getting Matrix Values...

## 2021-12-25 17:31:01 :

## 1 2 3 4 5 6 
## Plotting Embedding
## 1 2 3 4 5 6 7 8 9 
## ArchR logging successful to : ArchRLogs/ArchR-plotEmbedding-c6f5d16a29b-Date-2021-12-25_Time-17-31-01.log

# clean
p2MarkerHybridAll <- lapply(pMarkerHybridAll, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})
cowplot::plot_grid(plotlist=p2MarkerHybridAll)

Markers for hybrid1 cells

pMarkerHybrid1 <- plotEmbedding(
    ArchRProj = archProj, 
    colorBy = "GeneScoreMatrix", 
    name = markerHyrbid1[1:9], 
    embedding = "UMAP_cor",
    quantCut = c(0.01, 0.95),
    size=1/3,
    plotAs="points"
)

## Getting ImputeWeights

## No imputeWeights found, returning NULL

## ArchR logging to : ArchRLogs/ArchR-plotEmbedding-c6f51270ee1d-Date-2021-12-25_Time-17-31-09.log
## If there is an issue, please report to github with logFile!

## Getting UMAP Embedding

## ColorBy = GeneScoreMatrix

## Getting Matrix Values...

## 2021-12-25 17:31:10 :

## 1 2 3 4 5 6 
## Plotting Embedding
## 1 2 3 4 5 6 7 8 9 
## ArchR logging successful to : ArchRLogs/ArchR-plotEmbedding-c6f51270ee1d-Date-2021-12-25_Time-17-31-09.log

# clean
p2MarkerHybrid1 <- lapply(pMarkerHybrid1, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})
cowplot::plot_grid(plotlist=p2MarkerHybrid1)

Markers for hybrid2 cells

pMarkerHybrid2 <- plotEmbedding(
    ArchRProj = archProj, 
    colorBy = "GeneScoreMatrix", 
    name = markerHyrbid2[1:9], 
    embedding = "UMAP_cor",
    quantCut = c(0.01, 0.95),
    size=1/3,
    plotAs="points"
)

## Getting ImputeWeights

## No imputeWeights found, returning NULL

## ArchR logging to : ArchRLogs/ArchR-plotEmbedding-c6f550edb65c-Date-2021-12-25_Time-17-31-18.log
## If there is an issue, please report to github with logFile!

## Getting UMAP Embedding

## ColorBy = GeneScoreMatrix

## Getting Matrix Values...

## 2021-12-25 17:31:18 :

## 1 2 3 4 5 6 
## Plotting Embedding
## 1 2 3 4 5 6 7 8 9 
## ArchR logging successful to : ArchRLogs/ArchR-plotEmbedding-c6f550edb65c-Date-2021-12-25_Time-17-31-18.log

# clean
p2MarkerHybrid2 <- lapply(pMarkerHybrid2, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})
cowplot::plot_grid(plotlist=p2MarkerHybrid2)